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In this work we compute subleading oscillating terms in the Renyi entropy of Fermi gases and Fermi liquids 
corresponding to 2fcF-like oscillations. Our theoretical tools are the one dimensional formulation of Fermi 
liquid entanglement familiar from discussions of the logarithmic violation of the area law and quantum Monte 
Carlo calculations. The main result is a formula for the oscillating term for any region geometry and a spherical 
Fermi surface. We compare this term to numerical calculations of entanglement using the correlation function 
method and find excellent agreement. We also compare with quantum Monte Carlo data on interacting Fermi 
liquids where we also find excellent agreement up to moderate interaction strengths. 



I. INTRODUCTION 

Fermi liquids are an extremely common form of quantum 
matter that appear in a wide range of physical systems includ- 
ing alkali metals, overdoped cuprate superconductors, and 
high density quark matter (before color superconductivity sets 
in). Cousins of the Fermi liquid have been observed in the half 
filled Landau level HI and at a phenomenological level in lay- 
ered organic salts 1*21. The ubiquitous Fermi liquid has also 
played an important role in the recent exchange of ideas be- 
tween quantum many-body physics and quantum information 
science. This exchange has resulted in a deeper appreciation 
of the important role of entanglement, especially long-range 
entanglement, in the physics of quantum matter like Fermi 
liquids. Long-range entanglement is important because it dis- 
tinguishes interesting gapless phases like Fermi liquids and 
interesting gapped topological phases like fractional quantum 
Hall liquids from other more conventional phases like symme- 
try broken states that support long- range classical correlation 
but not long-range entanglement. 

Much of the discussion about entanglement has turned on 
a quantity known as the entanglement entropy. We consider 
a large quantum system AB divided into two components A 
and B (typically A and B are spatial regions, but other "en- 
tanglement cuts" are possible). The entanglement entropy is 
then the von Neumann entropy S{A) = — tr^(/9^ In {pa)) of 
the state of A where pA = Iib{pab)- When the state of the 
whole system, pab, is pure, the von Neumann entropy S{A) 
measures the amount of entanglement between A and B. We 
are especially interested in the case when pab — \G) {G\ with 
|G) the ground state of a local Hamiltonian. It has also been 
profitable to consider the Renyi entropy, 

5„(A)--^ln(tr(pl)), (1.1) 

1 — n 

which is actually a family of entropies labeled by n and which 
gives complete information about the spectrum of pA, the "en- 
tanglement spectrum". 

The basic rule governing the entanglement entropy in lo- 
cal ground states is the area law (for a review see Ref. ||3l). 
For a system with short-ranged interactions in d spatial di- 
mensions, the area law states that the entanglement entropy 



of a region A of linear size L grows like L'^^^, that is like 
the area \dA\ of the boundary dA of A. Fermi Uquids are ex- 
tremely interesting from an entanglement perspective because 
they possesses long-range entanglement that manifests as a 
violation of the area law |l4]-(8l. Indeed, entanglement entropy 
in a Fermi liquid ground state scales like L"^^^ In (L) hence 
showing a logarithmic violation of the area law. Gapless sys- 
tems in one dimension, including Luttinger liquids and quan- 
tum critical points, also show a logarithmic violation of the 
area law ||9l|T0l, while most other gapped and gapless systems 
in d > 1 dimensions obey the area law. In fact, the logarith- 
mic violations of the area law in d = 1 gapless theories and 
Fermi liquids are intimately related as we review below ||6l . 

Another important development in the study of many-body 
entanglement has been the appearance of numerical studies of 
entanglement in a wide variety of systems. Entanglement in 
one dimension has long been accessible using DMRG (for a 
review see Ref. flT\). and free fermions and bosons are acces- 
sible via the correlation matrix method |12|. More recently, 
quantum Monte Carlo and tensor network calculations have 
permitted computations of entanglement in simple quantum 
magnets |13 |, more complex spin liquid states|14|, topologi- 
cal states 1 15 1, and Fermi liquids |16|. We are thus finally in 
a position to begin a substantive comparison between theory 
and (numerical) experiment for universal terms in the entan- 
glement entropy. We now have agreement between theory and 
experiment in one dimension and for certain simple topolog- 
ical phases in two dimensions, and various other predictions 
have been validated at a more qualitative level, including a 
prediction of L In (L) entropy in a spin liquid with spinon 
Fermi surface f6', l4l and an observation of corner terms and 
terms associated with symmetry breaking in quantum mag- 
nets (although the agreement here is not yet precise) IfTTinrSl . 
Very recently the first calculations of Renyi entropy in inter- 
acting Fermi liquids were reported in Ref. llT6ll . with agree- 
ment at a quantitative level with previous theoretical predic- 
tions in Refs. |7 19] up to intermediate interaction strength. 
Although it should be noted that discrepancies between these 
calculations and the theoretical predictions increase in the low 
density limit of the Fermi liquid, this represents one of the first 
examples in more than one dimension of precise quantitative 
agreement between theoretical and numerical computations of 
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entanglement entropy in an interacting gapless system. 

Previous work has focused on the leading logarithmic term 
in the Fermi liquid which was argued to depend only on the 
geometry of the interacting Fermi surface and on the geom- 
etry of A. However, as may be expected on general grounds 
and as was evident in the data in Ref. lfT6l and elsewhere, 
there are also subleading oscillating terms in the Renyi en- 
tropy. Here we compute these subleading terms in the Renyi 
entropy analytically for the free Fermi gas. We also argue that 
the period of oscillation and exponent of the power law pref- 
actor are unmodified by interactions. We compare our results 
with extensive numerical data on free fermions and interact- 
ing Fermi liquids and find excellent agreement up to moder- 
ate interaction strengths. Thus we establish in considerable 
detail a quantitative agreement between theory and numerics 
regarding entanglement in Fermi liquids, both for the leading 
logarithmic violation as well as for the subleading oscillating 
term. 

This paper is organized as follows. We first review the one 
dimensional aspects of Fermi liquid entanglement before de- 
riving the form of the oscillating term in higher dimensions. 
We give an alternate proof for our formula when the region 
geometry is a long strip and further elucidate the structure of 
the entanglement spectrum in this case. Finally, we compare 
our theoretical predictions to numerical data for free and in- 
teracting Fermi Uquids and find excellent agreement. 



II. ENTANGLEMENT AND THE FERMI SURFACE 

We now describe the theoretical framework for our calcula- 
tions. Let R be the spatial region of interest in a Fermi liquid 
at zero temperature and let F denote the interacting Fermi sea. 
Once again, the Renyi entropy of region R is defined as 



Sn{R) 



1 



1 



ln(tr(p'A))- 



(2.1) 



It is well known that the leading behavior of the Renyi entropy 
is Sn ^ L"^^^ In L where d is the spatial dimension and L is 
the linear size of R IJHSl. Furthermore, the Widom formula. 



f \\ 1 

1 + - ^ 



provides a precise characterization of the prefactor of the loga- 
rithmic term in terms of the geometry of the interacting Fermi 
surface dV [5j-JJ. n^; and ri]^ are unit normals and the precise 
choice of length in the logarithm only modifies non-universal 
area law terms going like L'^~^. 

This formula can be obtained by describing the Fermi 
surface as an infinite collection of one dimensional modes 
(6|. Each such one dimensional mode is a gapless chiral 
fermion which contributes to the entanglement entropy like 
^-^^-^ where £ is some effective length (see below). Adding 
up these contributions for each point on the real space bound- 
ary dR and each mode on the Fermi surface dV leads to the 
Widom formula above. Even the dependence on the Renyi pa- 
rameter n is predicted by the theory since the Renyi entropy 



of a single interval in a CFT is flO]. Note also that 

each mode experiences many different effective lengths cor- 
responding to different one dimensional cuts through the real 
space region, but to logarithmic accuracy we may replace all 
such lengths with any particular representative of the linear 
dimension L. 

It is also known that the Renyi entropy has a subleading 
oscillating term in one dimensional Luttinger liquids. This 
term is analogous to Friedel oscillations and hence occurs at 
momentum 2kp 121]. For free fermions this term has the form 



Sn 



In 



COS {2kFl) 



with /3„ — 2/n and 



In 



r((i + n-i)/2) y 

l-n ^r((l-n-i)/2) ) 



(2.3) 



(2.4) 



(T{z) is here the gamma function) |21|. Below we will use 
our one dimensional formulation of higher dimensional Fermi 
liquid entanglement to demonstrate the existence of similar 
oscillating terms in the entanglement entropy in higher dimen- 
sions. 

First, however, we give a quick derivation of the Id result 
using conformal field theory. The Id electron operator may 
be decomposed as c{x) = CL{x)e^'^^''^ + CB{x)e'^^''^ where 
cl/h are the left and right moving halves a free relativistic 
fermion. They are the slowing varying fields entering the low 
energy description which we can use to compute universal 
terms in the entanglement entropy. The typical way we pro- 
ceed is to introduce twist fields or otherwise study the system 
on an n-sheeted surface to compute tr(p"). Here we simply 
note that an important feature of this setup is that translation 
symmetry is broken for n ^ 1. Furthermore, the singular 
branch points can in general produce localized relevant per- 
turbations bound to them. These perturbations produce cor- 
rections to the leading CFT scaling of Renyi entropy (there 
are also bulk irrelevant operators that contribute) since they 
change the free energy in the branched background. 

For example, the action is perturbed io S — Scft + 
g4>{zi) + g(j)(z2) where zi,2 are two branch points. The path 
integral is 



D(fields)e-'^[*''="'l 



ZcFTil - 9{{<t>{zi) + 0(2:2)))- 



-((</.(zi)+0(z2))' 



(2.5) 



Note that we have not been careful about the n dependence in 
this schematic expression. Since (0) = the first correction 
to the free energy (upon re-exponentiating the series and with 
the usual story about disconnected diagrams) is of order 
and is essentially a correlator of the form ((/)(zi)(/)(22))- To 
be specific, this is the part which depends on the separation 
of the branch points. We now go through this calculation in 
more detail for two branch points. 
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Consider a finite interval of length i in the free fermion 
CFT. To produce an n-branched surface with branch points at 
zi and Z2 with zi ~ Z2 = £, we can consider the conformal 
transformation given by 



z — Zl 
Z - Z2 



1/n 



(2.6) 



where w is a coordinate on a plane while z is the branched co- 
ordinate. Indeed, we see that winding z around zi or Z2 winds 
ui by a phase of 27r/n, and hence we must wind z by 27rn 
to wind w by 2tt. We now wish to compute correlations of 
operators inserted at zi and Z2 (the induced defect operators) 
to find corrections to the entropy. Given a primary field cf> of 
dimension A we wish to find 



(<^(z)<^(z')) = 



dw 
dz 



dw 
dz 



A 



(2.7) 



where the equality follows from conformal invariance under 
the transformation in Eq. |2.6| 

This equation is valid for all z and z' but we specifically 
want z = Zl + e and z' = Z2 + e with e a UV regulator 
Assuming ^ ^ e we obtain the following formulas; 



dw 
dz 

dw' 

17 " 



1/n 1 

e ' 

1/n 



and 



{(l){w)(t){w')) 



\w — w 



/|2A 



Putting everything together we find 

(0(zi)0(z2)) = ^ 



2A 



(D 



(D 



2A/r: 



2A/r: 



(2.8) 



(2.9) 



(2.10) 



(2.11) 



which explicitly shows a correction of the form ^-2A/n pj_ 
nally, to make contact with the Fermi gas result we must iden- 
tify the relevant operator, but this operator is just the 2kp den- 
sity operator given by e'^^'^'''^ c]^{x)cii{x) coming from the ex- 
pansion of c^{x)c{x). This operator has dimension A = 1 
and has an oscillating phase that is explicitly displayed. Thus 
the oscillating term term in the entanglement entropy is in- 
deed interpretable as a kind of 2k p density response, albeit in 
a branched space (which accounts for the strange scaling di- 
mension). Obtaining the prefactor for a given model requires 
more work since the field theory gives a cutoff dependent an- 
swer for the prefactor (which is anyway non-universal), but 
see Ref. |21 1 for a calculation for free fermions. 

We now return to the extension to higher dimensions. For 
concreteness, we set d = 2 and take i? to a disk of radius L 
and r to be a disk of radius kp. The extension of our results 
to arbitrary region shape (the analog of the Widom formula) 
is completely straightforward and will be recorded later First, 



we must be more careful about the effective length £ since 
we are studying subleading terms. The appropriate choice is 
£ ~ £{x, k) = 2L\ cos 9\, where 9 is the angle between and 
Uk, which is nothing but the chordal distance across the circle 
(parallel to Uk) at angle 9 \J_,^. This choice is the effective 
one dimensional distance experienced by a mode propagating 
in the direction starting at angle 9 on the circle and re- 
produces the correct thermal entropy to entanglement entropy 
crossover function Q|8]. 

The subleading oscillating term is then given by 



f J^l f /■ I cos {2kF£{x,k)) 

^"2^ 4 7,^ 4 {2kp£{xM^'- ■ 
Plugging in our various expressions we find 
AkpL 1 



(2.12) 



III 



8 {AkpL)^/" 



Tr/2 

X / d9{cos9y'^/" cosiAkpLcost 



(2.13) 



The final integral may be simplified when AkpL is large since 
the integrand is rapidly oscillating. Focusing on the part of the 
integral near 9 = Owe may write 



(2.14) 



where = ±7r/4 is an unimportant phase. Note that this 
leading estimate is completely independent of n as regards 
the integral over 9. Using this formula in our main expression 
gives, for AkpL large, the result 



(AkpL) 



1/2-2/ri 



COS {AkpL 



(2.15) 



where (j)' is another phase. Thus we have a prediction for the 
prefactor, power lay decay, and period of oscillation of the 
oscillating term in the Renyi entropy for a free Fermi gas. 



III. INTERACTIONS AND GENERAL REGION SHAPES 

When interactions are turned on we expect that the emer- 
gent C/(l)°° symmetry (the fact that the Fermi liquid fixed 
point is essentially free) will protect the exponent of the power 
law, but the numerical prefactor may be modified |7, 8|. 
For example, in one dimension the prefactor depends on the 
Luttinger parameter [21 1, but while the power law decay is 
also modified in one dimension, in higher dimensions the 
power law should remain unchanged since quasiparticles re- 
main sharp. The momentum of the oscillation will remain 
at the interacting 2kp, that is, while interactions may change 
the non-interacting Fermi surface, the correct momentum is 
always 2kp for the physical Fermi surface. 

The U{1)°^ symmetry refers to the fact that Landau's 
Hamiltonian for the Fermi liquid takes the form 



HpL = ^{^k- t^)nk + \'Y^ fkk' 



rikUk' 



(3.1) 



kk' 
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which has [HpL,nk] = for all k and hence an infinite num- 
ber of conserved quantities. Indeed, the ground state of this 
interacting Hamiltonian is a free fermion wavefunction with a 
Fermi sea F obtained by solving the self consistent equations 



and 



Efe < 0, fc e F 



e(fc) = Efc - ^ + ^ fkk' 
fc'Gr 



(3.2) 



(3.3) 



For example, if we work in a rotationally invariant system and 
impose Luttinger's theorem (by adjusting the chemical poten- 
tial) then the interacting Fermi surface always coincides with 
the free surface which is in turn determined solely by the den- 
sity. 

More generally, although the Fermi surface can change 
as interactions are tuned (e.g. in a lattice model with non- 
spherical Fermi surface), the entanglement entropy always 
obeys the Widom formula evaluated on the physical interact- 
ing Fermi surface. Similarly, all subleading corrections to the 
entropy in the toy model in Eq. 3.1 are those of a Fermi gas 
with Fermi surface F. However, we emphasize that this is 
not necessarily the correct answer for a physical Fermi liquid 
since the pure forward scattering model in Eq. |3.1| requires 
long-range interactions and hence is not in the same univer- 
sality class as a short-ranged Fermi liquid. This is reflected 

has an exact C/(l)°° 



3.1 



physically in the fact that while Eq. 
symmetry, the same symmetry in a short-ranged Fermi liquid 
is only emergent at low energies. 

Our argument for universality has recently been explicitly 
validated in Ref. Il22]| . That work examined a model of a 
Fermi liquid in which the quasiparticle residue could be con- 
tinuously tuned to zero at a quantum critical point, but where 
density and current correlation functions were those of a free 
Fermi gas. In the language of Fermi liquid theory, the Lan- 
dau parameters have been tuned to zero, but the quasiparticle 
residue is non-trivial. Nevertheless, the system is a non-trivial 
interacting Fermi liquid and Ref. |22| showed that the Widom 
formula for the leading logarithmic violation is still exactly 
obeyed. Furthermore, the oscillating terms we considered 
here are also present, unmodified, in the particular realization 
of an interacting Fermi liquid in Ref. Il22l . This suggests that 
the prefactor of the oscillating term can be expected to be near 
the free result provided the Landau parameters are small even 
if the quasiparticle residue is tiny. 

Now we turn to the general formula for the subleading os- 
cillating term in the Renyi entropy of a free Fermi gas for a 
convex but otherwise arbitrary region shape. Again, R is the 
real space region of linear size L and F is the spherical Fermi 
sea. We can also generalize to more complex Fermi seas in a 
straightforward way. Let x and k denote points on dR and dT 
respectively, and define the effective length £{x, k) to be the 
length of the line segment given by the intersection of the line 
{x + rtfcsjs S (— oo, oo)} and R. The convexity of R guar- 
antees that this intersection is a single line segment. It can 
also be checked that this definition reduces to our prescription 
for £ for a spherical region R given above. The subleading 



correction to the Renyi entropy is then 

A 1 / / I cos(2fcf£(a:,fc)) 



(3.4) 



We have already applied this formula to case when i? is a 
disk. It is also enlightening to consider a long strip region. 
Thus suppose the region i? is a long strip of length L and 
width W with L ^ W, and let 6 be the angle between the 
Fermi surface normal and real space normal. The effective 
distance is found to be ^ — M^/|cos6'| so long as ^ <C i 
which we take to be essentially infinite. Alternatively, we can 
consider a strip that wraps completely around the cycle of a 
torus of length L in which case translation invariance in the L 
direction is manifestly preserved. Within our approximation 
we find 



f^^kpL / cos( 

bn ^ :^ / dO COS 9 



27r 



ir/2 



\2kFW 



2/n 



COS 



2kFW 
COS 



(3.5) 

Assuming 2k pW ^ 1 we may perform the 6 integral as 
above by focusing on the region near = Q. The result is 



fnkpL ( 1 



27r \2kpW 



2/n 



2-K 



2kFW 



cos(2kFW). (3.6) 



Part of the reason why the strip case is interesting is that 
we may obtain the above result in another way. We can 
also say a great deal about the entanglement spectrum of the 
strip. The results are, however, restricted to free theories 
only. Thus consider again a free Fermi gas with spherical 
Fermi surface and examine the fermion two-point function 
G{x — y) = {c^ {x)c{y)) . This function may be obtained from 
a Fourier transform of the occupation number as 



G 



dPk 



ik-{x—y) 



(3.7) 



where = 9{kp — \k\). Our interest in this function is 
that, by virtue of Wick's theorem, it completely determines 
the reduced density matrix of a region R provided we restrict 
x,y ^ R. This is the standard correlation matrix method (see 
the appendix). 

Let K be the entanglement Hamiltonian (p = e^^) for the 
infinite strip (length L, see above). We know that K has the 
form 



^ ^ K{x,y)c'{x)c{y) 

x,y£R 



(3.8) 



because of Wick's theorem. Furthermore, the "matrix" K is 
related to G via G = -p^zp[ (proven by diagonalizing the "ma- 
trix" G with X, y restricted to R). Translation invariance in the 
L direction, call it the 2 direction, enables us to write 

K = K{x\y\x'',y^)^K{x\y\x^-y^O) (3.9) 

or upon Fourier transforming over — y^ 

E Kid{x\y\k)c\x\k)c{y\k). (3.10) 
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The new matrix Kid is related to the partial Fourier trans- 
form G. Consider the mixed position/momentum basis func- 
tion 



Gid(xi-yi,fc)=^G(a;i-yi,a;2)g- 



ikx 



27r 



e{kF - \/q^ + fc2)e*9(=^ -y ) (3.11) 



which is nothing but the two-point function of a one di- 
mensional Fermi gas with Fermi momentum k}^{k) = 
y/kj^ — fc2. Because both K and G are partially diagonalized 
by the Fourier transformation in the 2 direction we immedi- 
ately have that 



Gid(fc) 



1 



(3.12) 



Thus Kid is the single particle entanglement Hamiltonian for 
a one dimensional Fermi gas of Fermi momentum k]^ = 
y^kp ~ k'^ on an interval of length W. The full entanglement 
spectrum of the two dimensional strip is also now known in 
terms of the one dimensional spectrum e.g. we have a com- 
plete one dimensional spectrum for an interval of length W 
for each value of fc e [—kp,kp]. 

Using this information we can immediately check the lead- 
ing L log (W) term in the entropy. The entropy from each 
value of k is In (W) /3 and hence the total entropy is (to lead- 
ing order) 



S = L 



dk\og{W) kpL 



-kj 



2n 



In(W^). 



(3.13) 



A quick calculation with the Widom formula gives 

1 1 /■V2 



^=12 2^(2.L)(2M 



ir/2 



d6lcos6lln(Vl^) 



kpL 



\ii{W). 



(3.14) 



Returning now to the oscillating term, we can use Eq. 2.3 
to estimate the oscillating term for the strip in a different 
way. Each one dimensional spectrum identified above will 
contribute an oscillating term but with a variable k]^. The 
oscillating term is thus 



Sn ^ f n L 



dk C0s{2y^kj,- km) 

k^ 2n^2^/^^WwjP^' 



(3.15) 



As we have now repeatedly observed, if kpW is large then the 
integrand is rapidly oscillating and the integral is dominated 
by values of k near zero. Performing the effective Gaussian 
integral over k we find 



fnL nkp COS {2kFW) 



Sn 



fnkpL / 27r cos{2kFW) 



2tt V ^kpW (2fcj.M^)2/^ 



(3.16) 



(3.17) 



which is identical to our previous result. Note how the two 
methods obtain the same final form by integrating over either 
an effective length or an effective Fermi momentum. 



IV. COMPARISON BETWEEN THEORY AND NUMERICS 

We will now compare our theoretical predictions with nu- 
merical data for the Renyi entropy of the free Fermi gas for 
two 2D geometries, strips and circles. All numerical data 
was produced using the correlation function technique ex- 
tended to continuum Hamiltonians, and fits to the data were 
performed using a standard non-linear least squares (NLLS) 
package lfT2l [23] l24l . The data is presented in two plots with 
the corresponding parameters from the fits provided in the ta- 
bles. Computational details for computing Renyi entropies 
and the fitting procedure can be found in Appendix [A| 

For the circular geometry, we work in units where kp = 2 
so that the density of a spin polarized gas is n = kp/{4:7T) — 
1 /tt. This implies that the average number of particles in the 
real space circle is (N) = irL'^n = L^. In the figure we show 
Renyi entropies of a free Fermi gas computed numerically for 
disks with up to L = 6. The numerical data is fit to a two 
parameter function having the form 



Sn{L) = QiLln {L) + a2L + Qg ^2/1-1/2 ' ^4.1) 



with known parameters. 



ai 



(from the Widom fomula). 



n + 1 
in 



\fn\ 



27r 



81/2+2/n 



(4.2) 



(4.3) 



and 04 = 8. The parameters 02 and 05 are not predicted by 
our theory and are determined by fitting the numerical data. 
This form of the scaling law includes our new prediction for 
the subleading oscillations as well as the leading order term 
from the Widom conjecture and the area law term. In Table |l] 
we provide the parameters determined by our analysis as well 
as those from the fit. Root mean square residual (RMSR) val- 
ues, the square root of the sum of the residuals divided by the 
number of degrees of freedom in the fit, are provided as a met- 
ric for fit quality. Although we might expect other subleading 
terms not included in our theory to cause noticeable devia- 
tions, especially in the limit of small L, the fits remain close 
to the data over the entire range, as demonstrated in Figure |l] 
For the strip geometry we work with a slightly different 
density. We set our cell size to have length W = \/N. Work- 
ing in units where kj — 2^/tt, we fit to the form. 



Sn{W)^ai \n{W) + a2+az 



cos [a/iW + 05) 



^2/n+l/2 



Our analysis predicts 



ai 



n + 1 kfL 
in 27r 



(4.4) 



(4.5) 
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n 


ai 




aa 


a4 


(as) fit 


RMSR 




n 


Ol 


{a2)fit 


"3 


a4 


(as) fit 


RMSR 


1 


0.667 


1.774 











0.011 




1 


10.48 


32.95 











0.27 


2 


0.500 


1.384 


0.0253 


8 


0.60 


0.008 




2 


7.86 


25.58 


0.477 


7.09 


2.5 


0.24 


3 


0.444 


1.242 


0.0566 


8 


0.65 


0.015 




3 


6.99 


22.94 


1.026 


7.09 


2.3 


0.21 


4 


0.4167 


1.171 


0.0765 


8 


0.67 


0.027 




4 


6.55 


21.59 


1.357 


7.09 


2.3 


0.15 


5 


0.400 


1.123 


0.0869 


8 


0.69 


0.038 




5 


6.29 


20.76 


1.525 


7.09 


2.4 


0.26 



TABLE I. Comparison between theory and numerical data for the 
circular geometry. Parameters a2 and 05 are determined with a non- 
linear least squares fit while parameters 02,03,04, were fixed from 
theory predictions. The parameterization is given by equation |4.1| 
RA4SR is the root mean squared residuals from the non-linear least 
squares fit. The parameter 05 appears to be relatively insensitive of 
the Renyi order, while 02 does not. 
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FIG. 1. Renyi entropies Si through S^ for the 2D circular geometry. 
Solid lines are data and dashed lines are the fit whose parameters are 
found in Tablejl] This fit lines fall on top of the data for low order S„ 
and begin to deviate slightly for larger Sn ■ The amplitude and period 
of the oscillations, 03 and 04, are very accurately determined by the 
preceding analysis. 



(from the Widom fomula), 



03 



\fn\L^k 



2/n 



(4.6) 



and 04 = 2k f. We have chosen a fitting form for our data with 
two free parameters whose fits are provided in Table [III While 
the fits appear to be quite good there does appear to be small 
systematic differences between the theory and numerical data. 
This is due to a slight underestimation of the amplitude of the 
oscillations which may be an artifact of the approximation to 
the integral over ( 



' in equation 



3.6 and 3.15 



However, the am- 
plitude is qualitatively correct and the period is accurate. We 
expect that we can improve agreement between our analysis 
and numerics by increasing W , which is currently limited by 
computer memory. We note that direct comparisons between 
the quality of fit metrics for the two geometries and different 
order Renyi entropies need to be considered carefully as the 
scales of the entropies are not the same and the RMSR values 



TABLE II. Comparison between theory and numerical data for the 
strip geometry. Parameters 02 and 05 are determined with a NLLS fit 
while parameters 02,03,04, were fixed from theory predictions. The 
parameterization is given by equation |4.4| Similar to the circular 
geometry, parameter 05 appears to be independent of Renyi order 
while 02 does not. 
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FIG. 2. Renyi entropies Si through Ss for the 2D strip geometry. 
Solid lines are data and dashed lines are the fit parameterized in Table 
[n] Here we note good agreement with the period, 04, and acceptable 
agreement for the amplitude, 03, with our analysis. 



are scale dependent. 

Finally, let us mention what happens with interactions. For 
this we may consider the data in Ref. fl^. In that work Renyi 
entropies were computed for various interacting Fermi liquids 
as a function of interaction strength for both short-range in- 
teractions and Coulomb interaction using the swap operator 
in quantum Monte Carlo|13|. For 6*2, even at moderate in- 
teraction strengths (up to = 1 with Coulomb interactions), 
the free and interacting Renyi entropies agree quite well. This 
data thus provides numerical evidence for our earlier theoreti- 
cal claim that the period and power law prefactor of the decay 
is unchanged by weak to moderate interactions. Additionally, 
as we described above, Ref. 1221 has shown that for a certain 
solvable model of an interacting Fermi liquid, the proposed 
universality of Renyi entropies is indeed exactly true. 



V. CONCLUSIONS 

In this paper we have considered subleading oscillating 
terms in the Renyi entropy of free Fermi gases and Fermi 



7 



liquids. We give a simple analytic formula for the sublead- 
ing oscillating term using the one-dimensional formulation of 
Fermi liquid entanglement. This formula compares favorably 
with numerical calculations of the Renyi entropy of free Fermi 
gases in various geometries. We also considered the effects of 
interactions and argued for a certain degree of universality in 
the subleading oscillating term. Our arguments were checked 
by comparing to previous quantum Monte Carlo calculations 
of Renyi entropies in Fermi liquids as well as by comparing to 
exact results in a solvable model of a Fermi liquid. Thus we 
have established precise agreement between theory and nu- 
merics for both the leading and subleading terms in the Renyi 
entropy of Fermi liquids. 

An interesting future direction, which we have only just 
touched on here, is the exploration of the physics of the entan- 
glement spectrum. We know the full spectrum for a strip, but 
more generally it would be desirable to have an understand- 
ing the spectrum of more general regions and in the presence 
of interactions. Something like the bulk-edge correspondence 
for topological phases should be valid for Fermi liquids as 
well, but the precise form of this correspondence remains un- 



certain. It would also be very interesting to extend our results 
to other kinds of quantum matter which support a Fermi sur- 
face but which may not be simple Fermi liquids. Entangle- 
ment entropy in these models provide a conceptually simple 
probe of the Fermi surface, even if it is not associated with 
conventional electrons, and since both the oscillating term and 
the leading logarithmic term know about the Fermi surface 
geometry, we can extract the universal prefactor (analogous 
to the central charge) in front of the logarithmic term in the 
Renyi entropy. Finally, although we considered only spheri- 
cal Fermi surfaces here, it is possible to extend our results to 
more general Fermi surface shapes. 
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Appendix A: Computational Details 

We compute the Renyi entropies from the eigenvalues of a 
spatially reduced density matrix using the correlation function 
technique Iil2i . For a free particle Hamiltonian, 

H = — ^ ^ im,n(^n^m (Al) 

with n and m subsystem site indices, and ci and c| the creation 
and annihilation operators for state i, eigenvalues for the den- 
sity matrix can be computed using the relationship between it 
and the correlation function matrix. The correlation function 
matrix, 

= (cjc,), (A2) 

is determined entirely by the one particle operators. For this 
Hamiltonian we can write the density matrix as the exponen- 
tial of a fictitious Hamiltonian, the entanglement Hamiltonian, 
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with energy levels and single particle operators ak and a^, 



p = K. exp - ^ alak^k 



(A3) 



where /C is a normalization constant set by Tr(p) = 1. The 
new states, are related to the eigenvectors of the correla- 
tion matrix by, c; = J2k 4'k{i)'^k, and the eigenvalues of the 
entanglement Hamiltonian are related to those of the correla- 
tion matrix, Q, by. 



C« = log 



(A4) 



The von Neumann Entropy, ^i, is then computed as. 



5i-^ln(l + exp(eO) + J 



-6 



+ cxp {-^i) 



(A5) 



and higher order Renyi entropies using a recursive formula. 



Too = 

Wo = l 

i = 1 

while i < Ntot 

Wi = Wi-i{l + exp(-^j)) 

TOj = TOi_i + exp + TOi_i) 

i = i + l 



(A6) 



as shown in 1*231 • 

When computing high order Renyi entropies or the entropy 
of a large region, to and w can develop numerical instabilities. 
The eigenvalues of the correlation function matrix, Ci, can be 
interpreted as the occupations of orbitals for region A [251 . 
As the region size increases more of these orbitals become 
occupied with Q 1 exponentially fast. The normalization 
factor, w, and the unnormalized trace, m, contain terms that 
diverge as the system or Renyi entropy order, n, grows large: 

w = a(i - Q)-' and TO ^ n.(i - c.)-"- 

While each term diverges individually, the ratio, r„j^ = 



("'"to 



0. This ratio can be written as a function of the 
Renyi entropies which has a known scaling form, r„ = 
exp ((1 — n)Sn)- We can see that this ratio, rm.w, scales as 
W^^" for cylindrical geometry and l'^^")^ for the spher- 
ical geometry, both of which go to zero as the system size 
increases. To deal with these numerical issues and maintain 
good accuracy for all system sizes, we use arbitrary precision 
libraries. 

For the fits, very small L and W data are excluded and re- 
gion sizes spanning less than one half the box size were used. 
This eliminates the contributions from region sizes which 
have very few particles on average and minimizes the interac- 
tion between periodic images of the subsystem. We fixed the 
parameters given by our analysis and fit the unknown param- 
eters. The fitting procedure was performed using a non-linear 
least squares fit which minimized the sum of the residuals. We 
note that fits were also made without fixing the known param- 
eters to their theoretical values, however the resulting values 
were determined to be too sensitive to the range of data in- 
cluded in the fit. 



